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Abstract 



We consider an array of fiber lasers coupled through the nearest neigh- 
bors. The model is a generalized nonlinear Schroedinger equation where 
the usual Laplacian is replaced by the graph Laplacian. For a graph with 
no symmetries, we show that there is no resonant transfer of energy be- 
tween the different eigenmodes. We illustrate this and confirm our result 
on a simple graph. This shows that arrays of fiber ring lasers can be made 
^T J temporally coherent. 



1 Introduction 



The dynamics of coupled oscillators continues to be an active area of research 
due in large measure to its ubiquitous presence in a wide range of disciplines 
such as social, life and neuro-sciences [T] . It also applies in the traditional areas 
of classical and quantum mechanics, in particular in nonlinear optical systems 
describing phenomena such as light localization in waveguide arrays. 

This work concentrates on a simple model to understand how synchroniza- 
tion can be achieved when light is propagating in ideally identical optical fibers 
coupled by a suitable scheme. Two well known models encapsule most of the 
behavior considered here: (i) The Kuramoto (KM) model 

— * = m + K}^ sin(4>j - 4>j) (1) 
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(ii) The discrete nonlinear Schrocdinger (DNLS) model 

i-Tf- = c(u l+ i + +Ui-i) + \ui\ 2 Ui (2) 

where the latter is the most relevant to our study. The collective (synchronous) 
behavior of these coupled optical fibers would produce a maximum power elec- 
tromagnetic output. Historically, research in optical fiber technology has been 
driven by the increased demand for better communication systems and a natural 
by-product has been the development of fiber lasers 3 ; . These are interesting for 
applications like for example material cutting where they are starting to replace 
solid state devices. To increase their power output, several groups are consider- 
ing arrays of fiber lasers. Here light is amplified in individual, decoupled fiber 
amplifiers which are then coupled within a single cavity or in a second cavity. 
The grand challenge is to achieve a coherent power output that scales with the 
square of number of elements in the array. In this context, temporal coherence 
is the equivalent of synchronization of oscillators in an array. In practice, this 
synchronization has proven to be very difficult. Typically the efficiency dra- 
matically diminishes as the number of fiber amplifiers increases. Amongst the 
various combining schemes a most interesting case is that of Fridman et.al [4]. 
There the authors combine passively as many as 25 fiber lasers in a two dimen- 
sional array. They analyze the power output in a short time interval over many 
round trips and observe two main effects. First, even though the efficiency of 
the phase locking is around 20 to 30%, there are rare events where it exceeds 
90%. The second observation is that the efficiency depends on the coupling 
architecture, the best performance is when connectivity is increased. 

These fiber arrays at high intensities are well described by the DNLS equa- 
tion @. A number of studies of this model suggest that nonlinearity can 
enhance coherence. Following this, here we propose a more general coupling 
scheme where the usual discrete Laplacian operator in the DNLS equation is 
replaced by a graph Laplacian. Our main result is that when choosing a specific 
form of the coupling we are able to "separate" the modes so the nonlinearity 
will only couple them weakly. Then the system can be considered as temporally 
coherent. This holds for any graph such that the eigenvalues of the Laplacian 
are distinct, i.e. a graph with no symmetries. We illustrate this by analyzing 
a simple graph and confirm this result numerically. We calculate the coherence 
factor of the array for different initial mode distribution and show that fiber 
ring lasers in an array can be synchronized. 

The article is organized as follows. We introduce the model in section 2 and 
show that the nonlinearity does not couple the linear eigenmodes on average so 
that their average amplitude is constant. In section 3, we illustrate this general 
result on a particular graph. This system is analyzed numerically in section 
4. There we confirm our predictions and consider different mode distributions. 
Conclusions are presented in section 5. 
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Figure 1: The top panels presents two physical arrangements of a four fiber 
array. On the left the fibers 2 and 3 are closer and thus coupled stronger than 
the other pairs of fibers. On the top right panel the couple 3 and 4 are closer. 
The bottom panel shows the graph used to model these devices. The arrows 
on the branches are oriented arbitrarily and the coupling parameters a, /3 are 
adjusted to represent each of the physical configurations. 

2 The model 

We propose a design of a fiber array where the coupling is purposely made 
irregular. We expect that in a network the coupling will be different due to 
the arrangement of the single fibers in relation to others. Fig. [1] presents the 
cross-section of the distribution of fibers in the array and the associated network 
or graph used to model it. A consequence of this is that we need to generalize 
the classical 2D nonlinear Schroedinger equation 

iut = Au + \u\ 2 u , 

to the graph NLS 

iu t = Gu + \u\ 2 u, (3) 

where the standard 2D Laplacian is replaced by the graph Laplacian G. This 
operator is a natural extension of the continuous Laplacian see [5] for examples 
of its application. The specific form of ([3]) written for the graph shown in Fig. 
[1] gives 

iu\ = (u 2 - «i) + |ui| 2 wi , 

iu% = - (u 2 - Ui) + a (u 3 - u 2 ) + {ua - u 2 ) + \u 2 \ 2 u 2 , 

iii 3 = ~a (u 3 -u 2 ) + f3 (u 3 - u 4 ) + |u3| 2 M3 > 

iu 4 = (u 2 — u 4 ) — P (u 3 — ua) + |u4| 2 «4, 



where Ui, i = 1 — 4 are the values of the field amplitude, where ft is the coupling 
for the branch between nodes 1 and 2, where a is the coupling for the branch 
between nodes 2 and 3 and where 1 is the coupling for the branch between nodes 
2 and 4. The system of equations above can be written symbolically as 



iU t = GU + N(U). 



(4) 



In the above equation the vector U , the graph Laplacian G and the nonlinearity 
are respectively 
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This system preserves the "mass" 



M =T 



Uj 



and the energy 



h = J2W+J2 



N 4 , 



where V is the discrete gradient associated to the graph [5]. 

Since the matrix G is symmetric it is natural, following [5], to use as a basis 
for U the eigenvectors of G, such that Gzi = —ufzi. These are the columns of 
the orthogonal Z matrix such that 



GZ = ZD, 



(5) 



where D is the diagonal matrix of diagonal —uj 2 , —uj 2 ---. We then introduce the 
vector r of components 7^, k = 1, . . . n such that 

U = ZT. 

In terms of these new coordinates jk , equation (j4]) reduces to 



Hk 



-Wfc7fc 






(6) 



where we have used the orthogonality of the eigenvectors of G, Z l = Z T . The 

E* 



term jUjpMj can be written as 



\Uj\ 2 Uj 



We then get the final equation 

iik 



ullk 



■ 2J ZjlZj m Zj n JlJ m J^ 
jlmn 



(7) 



This is the equation that we will analyze throughout the article. 

Let us assume that all the eigenvalues ujj are simple. This is the generic 
case for a graph without symmetries. Then we can simplify the equation (|7|) by 
eliminating the first term on the right hand side. For that we introduce 

Ik = e^%, 
to obtain 

jlmn 

Only the resonant terms such that 



u 



2 , ,2 



loI = 0, (9) 



will contribute on the long term to /3k- Then, because the frequencies are all 
different, the resonant condition is satisfied if I = n and m = k. We then obtain 
the resonant evolution of (3k as 



H3k = (3kY z Ji z i*\h\ ■ ( 10 ) 

ji 



This equation is such that the intensity in each mode is constant 

h = \(3 k \ 2 , 4 = 0. 
The solution of equation ([TO)) is then 

Pk(t) = p k (0)exp{-i t n k ), (ff) 

where the nonlinear correction to the frequency is 

n k = Y / 4 z ^\Pi\ 2 - (12) 

Returning back to the original variables jk we get our main result 

7fc(*) = 7fc(0)exp [i t (ojI - n fc )] . (13) 

There is no long term energy transfer between the modes. For a given initial 
condition, the modal distribution is fixed and the system will keep this for all 
times. 

Before presenting a case study, a few remarks are in place. Our analysis is 
valid for any graph, of arbitrarily large size. For our approximate analysis to be 
valid it is important that eigenvalues be different and well separated. Then the 
phase resonance condition will be well satisfied. We will illustrate this below 
by comparing our prediction to the numerical solution of the full problem. Also 
note that when there are multiple eigenvalues, the change of variable from 7 to 
(3 is no longer valid and the degeneracy will induce a chaotic temporal dynamics. 
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Table 1: Eigenmodes (Xi,Zi) for the tree. 

3 A case study: a graph with a swivel 

We now illustrate the above considerations on a simple example derived from the 
network shown in Fig. [TJ This will enable us to put numbers on the frequencies 
uj k and f2 fc . In the study [S], we introduced the notion of "soft node" where 
the eigenvector has a zero coordinate. These are important because if applied 
to them the damping or forcing of the network is ineffective. We will consider 
the simplest such network that has a soft node. It corresponds to the tree with 
(3 = 0. 

The eigenvalues and eigenvectors Zi (columns of the z matrix) can be com- 
puted analytically, they are given in Table 1 in terms of a, 



5 = \/4a 2 - 4a + 9, 



and 



n 3 = y/2 + (6-2a- l) 2 /4 + (6 - 2a + 3) 2 /4, 



n 4 = ^2 + {S + 2a+ l) 2 /4 + (6 + 2a - 3) 2 /4. 
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Figure 2: Plot of the eigenfrequencies Ui, i — 1 — 4 as a function of a for the 
graph shown in Fig. [1] 
with /3 = 1. 

The mode Z2 which is independent of a is shown schematically in Fig. [3] 
where we have plotted the magnitude and sign of the coordinate using vertical 
arrows. The arrows are opposite and equal for z\ and z| because the nodes 1 
and 4 play a symmetric role i.e. the graph is invariant by the automorphism 
transforming node 1 to node 4 [6]. 




Figure 3: Schematic representation of the constant mode z-z-, corresponding to 
U)2 = 1 
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Figure 4: Plot of the components of the non trivial normalized eigenvectors z 3 
(left panel) and z 4 (right panel) as a function of a. 

Let us compute the frequencies of oscillations 17 & of the resonant modes. The 
right hand side of relation (jllj) can be written in matrix form 



mi = n 



where the vector I = 



h 
h 

w 



with Ij = \(3j\ 2 , where the vector 57 = 



where the matrix M has a general term 



\n 4 J 



and 



m k i 



E 



Zjk Zjl- 



(14) 



Introducing the matrix Z from Table 1 and computing M using Matlab we get 
the general linear system 

/0.5 0.5 0.5 0.5 \ //A /fiA 

0000 i 2 _ n 2 

m 32 m 33 m 34 J 3 3 

y m 42 m 43 m 44 y \/ 4 / \fi 4 / 

This matrix enables to compute the nonlinear corrections Qk once the mode 
distribution 7^ is given. The matrix M has rank 3 and its image is 

Im(M) = {x,z,t\. 

Note that 51 2 = 0. Taking for example I\ — 2, 7 2 = ^3 = I 4 = we get 51 x = 
and 57 2 = 51 3 = 51 4 = 0. 



4 Numerical results 

To validate our approach we now solve numerically the original equations Q for 
the tree configuration described above ((3 — 0). We chose the initial mode con- 
figuration and let it evolve. The resolution was done using the ode45 subroutine 
of Matlab [7]. 



Fig. [5] shows the time evolution of the energies Ik for k = 1 — 4 for two 
values of the coupling, a = 0.5 (left panel) and a — 3.5 (right panel). The 
initial conditions are the same 0i = 0,(32 = 1-3, /3 3 = and P4 = 2. Notice 
how I4 is approximately constant in agreement with our prediction. On the 
other hand in both cases the modes 2 and 3 have close frequencies so they 
couple strongly. For a = 0.5 the frequencies are closer so that more interaction 
occurs. For a = 3.5 we have W4 >> u>2 so the phase term in ([9]) is larger and 
the approximation is better. Furthermore, it is clear that for this case a slow 
dynamics will show a fixed point for mode 4; and homoclinic orbits for modes 
1-3. All this is seen in the right panel of Fig. [5j 
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Figure 5: Time evolution of the energies in each mode Ik for k = 1 — 4 for two 
values of the coupling, a = 0.5 (left panel) and a — 3.5 (right panel). 

To relate our results to figures of merit in experimental work on coherent 
beam combining, and in particular paying attention to that of Fridman et al. 
[4], we introduce the coherence factor C 



C 



E, 



E, 



(15) 



For the modes w,:, i — 2— 4 this coherence factor will be smaller than one because 
the eigenvectors Zi have to have positive and negative coordinates, which could 
be interpreted as being in (same sign) or out (different sign) of phase data. As 
an example consider that 



7i = 72 = 73 = 0, 74 ^ 



Then for a large we have 



"2 = ^2474 



-O.874 



W 3 



-23474 ~ O.674 



and u\ = U4 = 0. This clearly represents an out of phase combination leading 
to a low coherence coefficient value, 



c= u 2 + u, 1 



U2 



U4, 



If however one considers instead of the nodes u, the modes 7, say by applying 
the linear transformation given by the matrix Z, then the coherence coefficient 
becomes 1. This could be done using transformational optics. 

Another direction would be to use the Goldstone mode z\ corresponding to 
all node components m equal. For this we need to guarantee that the mixing 
with the other frequencies remains small. This can be done for the tree as we 
now show. We fix 71 = 2, 74 = 0, and choose two sets of the pair (72,73), 
72 = 0.1, 73 = and 72 =0, 73 = 0.1 . The results are shown in Fig. [5] with 
72 = 0.1, 73 = on the left panel and 73 = 0.1, 72 = on the right panel. As 
expected on the left panel we see a stronger coupling than on the right panel 

because 

2222 
uj 1 -u) 2 < uj 1 -u> 3 , 

so that the rotating wave approximation leading to (|13[) is better satisfied. Con- 
sidering the coherence factor, the parameters 71 = 2, 74 = 0, and 73 = 0.1, 72 = 
will give CkI. 
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Figure 6: Time evolution of the energies in each mode Ik for k = 1 — 4 for 
71 = 2, 74 = 0, . The left panel is for 72 = 0.1, 73 = and the right panel for 
73 = 0.1, 72 = 0. The other parameter is a = 3.5 

The results of this section show that if the array is "prepared" in a given 
mode, it will stay on that mode on average. For this preparation, one can force 
the array at resonance like in [5]. As an example, for our tree configuration we 
select mode 4 by damping nodes 1 or 4 and forcing only nodes 2 or 3. 

5 Conclusions 

Motivated by studies in fiber laser arrays and how the different coupling schemes 
between fibers can affect the coherence property of the output, we introduced 
an array where the coupling between fibers is irregular. This leads to a discrete 
nonlinear Schroedinger equation where the usual discrete Laplacian is replaced 
by a graph Laplacian. 

To connect our work with experiments, we computed a coherence coefficient 
and showed that we can synchronize the array. This study then opens the way 
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of using a new approach based in graph Laplacians for designing new coupling 
schemes for arrays of fiber lasers. 
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